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On the basis of the Bogohubov-de Gennes theory we study the transformation of the quasi- 
particle spectrum in the mixed state of a mesoscopic superconductor, governed by an external 
magnetic field. We analyze the low energy part of the excitation spectrum and investigate the 
field-dependent behavior of anomalous spectral branches crossing the Fermi level. Generalizing the 
Caroli-de Gennes-Matricon approach we present an analytical solution describing the anomalous 
branches in a vortex with an arbitrary winding number. We also study the spectrum transformation 
caused by the splitting of a multiquantum vortex into a set of well separated vortices focusing mainly 
on a generic example of a two-vortex system. For vortices positioned rather close to the sample 
surface we investigate the effect of the quasiparticle reflection at the boundary on the spectrum and 
the density of states at the Fermi level. Considering an arbitrary surface curvature we study the 
disappearance of an anomalous spectral branch for a vortex leaving the sample. The changes in the 
vortex configuration and resulting transformation of the anomalous branches are shown to affect 
strongly the density of states and the heat conductance along the magnetic field direction. 

PACS numbers: 74.25.0p, 74.78.Na, 74.25.Fy 



I. INTRODUCTION 

Modern technology development provides a unique 
possibility to study exotic vortex states in mesoscopic 
superconducting samples of the size of several coherence 
lengthsii^i^. Tuning the external applied magnetic field 
one can switch between a rich variety of energetically 
favorable or mctastable vortex configurations which can 
not be realized in bulk systems. Of particular interest is 
a possibility to obtain multiquanta (giant) vortex states 
with winding numbers larger than unity for certain in- 
tervals of external magnetic field (see, e.g., Ref. |2). The 
merging of individual vortices into a multiquantum one 
occurs under the influence of screening currents which 
push the vortices to the sample center. Note that al- 
ternatively the stable multiquanta vortices can appear 
even in a bulk superconductor because of the pinning on 
columnar defects with radii of the order of the coher- 
ence length^. Experimentally the vortex configurations 
in mesoscopic systems and the phase transitions between 
them can be studied, e.g., by Hall probe measurements 
of the branching of the magnetization curved or by ob- 
servation of the vortex entry into the sample using the 
point contact techniques^. 

Changing the number and arrangement of the flux lines 
we can tune the low energy excitation spectrum which is 
known to be responsible for low temperature thermody- 
namic and transport properties of the sample. The mech- 
anism of such changes in the subgap quasiparticle spec- 
trum is associated with the modiflcation of the anoma- 
lous energy branches crossing the Fermi level. For well 
separated vortices positioned at distances much larger 
than the core radius the behavior of the anomalous bran- 
ches can be described by the Caroli-de Gennes-Matricon 
(CdGM) theory^. For each individual vortex the en- 
ergy e(/i) of a subgap state varies from — Aq to +Ao as 



one changes the angular momentum /i defined with re- 
spect to the vortex axis. At small energies |e| ^ Aq the 
spectrum is a linear function of ^: e(^) — — /ii^, where 
Lo Si Ao/(A:j^^), Aq is the superconducting gap value far 
from the vortex axis, k±_ = \Jk\ — k\, kp is the Fermi 
momentum, fc^ is the momentum projection on the vor- 
tex axis, ^ = fi-Vp/Ao is the coherence length, Vp is the 
Fermi velocity, and jj, is half an odd integer. 

With the decrease in the intervortex distance the 
quasiparticle tunneling between the vortex cores comes 
into play resulting in the modification of the anomalous 
branches^. Finally, when the vortex cores merge one 
obtains a multiquantum vortex with a certain winding 
number M. The number of anomalous branches per spin 
projection^ is conserved during this process of crossover 
from M individual flux lines to the M-quantized giant 
vortex and equals to the vorticity M. Previously, the 
behavior of the anomalous branches in a multiquantum 
vortex has been investigated numerically^ and analyti- 
cally for a step-like model profile of the order parameter 
in the core^^. For vortices with an even vorticity all the 
anomalous branches cross the Fermi level at nonzero im- 
pact parameters b ~ —i^i/kj_: 



e{^i)^-{pi±^i,)Ao/{k^O 



(1) 



where j = I...A4/2, fiM/2 ~ k±^. For a vortex with an 
odd winding number there appears a branch crossing the 
Fermi level at zero impact parameter. 

The wave functions of the subgap states are localized 
inside the vortex core because of the Andrcev reflection 
of quasiparticles at the core boundary. Any additional 
normal scattering process should modify the behavior of 
the anomalous spectral branch. Such modification can be 
caused even by atomic size impurities, as it was predicted 
by Larkin and Ovchinnikov in Ref. I 111 . For a vortex ap- 
proaching a flat sample boundary the distortion of the 



local density of states (DOS) profile has been analyzed 
in Ref. [ij numerically on the basis of Eilenberger theory 
both for s- and d- wave pairing symmetries. Certainly 
the role of normal scattering at the boundaries can be of 
particular importance for vortices trapped in mesoscopic 
samples. For a single vortex placed in a superconduct- 
ing cylinder an appro pria te spectrum transformation has 
been studied in Refs. [iSlflJI . 

Experimentally the behavior of the anomalous bran- 
ches can be probed, e.g., by the scanning tunneling 
microscopy (STM) or by the heat transport measure- 
ments. The modern STM techniques is a unique tool 
for the study of the local DOS profiles and, thus, could 
provide us the information about the number and con- 
figuration of the spectral branches crossing the Fermi 
level. An important advantage of the heat conductance 
measurements along the vortex lines is associated with 
the fact that probing the number and transparency of 
quasiparticlc transport channels this method appears to 
be sensitive also to the fc^-dispersion of the spectrum 
and, in particular, to the group velocity of quasiparti- 
clc modes propagating along the vortex core o^^i^^ . In- 
deed, it is the small group velocity of CdGM states which 
is responsible for a strong suppression of the heat con- 
ductance Kv '■^ T^kp^/ifiAo) along a singly-quantized 
vortex core as compared to the Sharvin's conductance 
^S/i ~ Tikp^Y' /h of a normal metal wire of the radius ^ 
at a certain temperature T: 



Ky IT 



<1 . 



(2) 



Within the Landauer approach such suppression of the 
vortex heat conductance can be understood as a conse- 
quence of a strong reduction of the effective number of 
conducting modes N^ = k^/kq-, where kq = 7rT/(3?i) 
is the universal heat conductance per conducting mode 
in a normal metal. Taking the interlevel spacing too ~ 
Ao/(fci?0 for the anomalous branch at k^ = 0, one ob- 
tains Ny ~ T/ujQ which agrees with the above estima- 
tion ^. Both the intervortex quasiparticlc tunneling 
and the normal scattering at the sample boundary are 
expected to affect the effective number of conducting 
modes resulting in the dependence of the heat transport 
on the vortex configuration in a mesoscopic superconduc- 
tor. The increase in the heat conductance stimulated by 
the boundary effects was demonstrated in Ref. [ij for a 
single-vortex state in a cylinder. 

It is the goal of the present paper to study both the 
transformation of the anomalous branches and distinctive 
features of the DOS and heat transport in different multi- 
vortex configurations in mesoscopic samples. We include 
in our consideration both the spectrum transformation 
caused by the giant vortex splitting and the process of 
an anomalous branch formation (disappearance) which 
occurs when a vortex enters (exits) the sample. In the 
present study we address only the case of homogeneous 
mesoscopic superconductors without any defects or pin- 
ning centers. 



The paper is organized as follows. To elucidate our 
main results we start from a qualitative discussion of the 
behavior of the anomalous branches in a mesoscopic su- 
perconductor (see Sec. |ll]). In Sec. IIIII we introduce 
the basic equations used for the spectrum calculation. In 
Sec. IIVI we consider an analytical solution describing the 
spectrum of a multiquantum vortex line. In Sec. |V] we 
study a generic example of the spectrum transformation 
caused by the decay of giant vortices, i.e. the splitting 
of a doubly-quantized vortex into two individual singly- 
quantized vortices. The influence of the normal reflection 
at the sample surface on the quasiparticlc spectrum for 
a vortex positioned close to the boundary is analyzed in 
Sec. Ell In Sec. Ellland Sec. Ivml we calculate the den- 
sity of states and the thermal conductance, respectively, 
using the quasiparticlc spectra found in the previous sec- 
tions. We summarize our results in Section [IXI Some of 
the details of our calculations are given in appendices. 



II. TRANSFORMATION OF ANOMALOUS 

SPECTRAL BRANCHES: QUALITATIVE 

PHYSICAL PICTURE 

As we discussed in Introduction there are two basic 
mechanisms responsible for the transformation of anoma- 
lous branches in a multivortex configuration: (i) the tun- 
neling of quasiparticles between the vortex cores; (ii) 
quasiparticlc scattering at the sample boundaries which 
comes into play when the vortices approach the super- 
conductor surface. To clarify the key physical conse- 
quences of these mechanisms hereafter we consider two 
model problems: (i) electronic structure of a multivor- 
tex system positioned rather far from the boundary; (ii) 
electronic structure of an individual vortex approaching 
the sample surface. 



A. Effect of intervortex quasiparticle tunneling 

Let us start with a qualitative analysis of the behavior 
of the anomalous branches and consider a set of vortex 
lines parallel to the z-axis. It is the case of intermedi- 
ate values of magnetic field when vortices are quite far 
from the boundary but do not merge into a multiquan- 
tum vortex. In the (a;y)-plane the vortex centers defined 
as points of the order parameter phase singularities (and 
hence as zeros of the superconducting order parameter) 
are positioned at certain coordinates r^. 

The system is homogeneous in the z-direction and, 
as a result, the momentum kz appears to be conserved. 
The two-dimensional quantum mechanical problem in 
the (xy)-plane can be strongly simplified provided the 
wavelength fc]] is much less than the superconducting 
coherence length ^. Thus, following standard quasiclas- 
sical procedure (see Section IIIII for details) one can de- 
scribe the quantum mechanics of quasiparticles using the 
geometrical optics picture. An important distinctive fea- 



turc of this picture in superconductors is that all the 
classical trajectories can be approximately considered as 
straight lines. The bending of these straight trajecto- 
ries is negligible because of a small momentum change 
Sk ^ 1/^ during the process of quasiparticle scattering 
at the inhomogeneous superconducting gap profile, i.e. 
during the Andrccv reflection. We also can neglect the 
trajectory bending caused by magnetic field, since we as- 
sume the cyclotron radius rn ^ Vf/ujh ^ kp£,'^{Hc2/H) 
to exceed all the relevant length scales of our problem. 
Here ujh = \e\H/(mc) is the cyclotron frequency, m is 
the electron effective mass and Hc2 is the upper critical 
field. 

For quasiparticlcs propagating along the classical tra- 
jectories parallel to kj^ = k±(cos 9p, sin 9p) we introduce 
the angular momenta fi = [r, kj^] • Zq = k±r sm{9p — 9) 
and jli ~ fi — [vi, \ij_] ■ zq defined with respect to the z- 
axis passing through the origin and with respect to the 
i-th vortex axis passing through the point r^ , correspond- 
ingly [{r,9,z) is the cylindrical coordinate system]. Ne- 
glecting the quasiparticle tunneling between the vortex 
cores and the normal scattering at the sample boundary 
we get degenerate CdGM energy branches: Si = —uojli 
for |e| ^ Aq. For a fixed energy e we can define a 
set of crossing quasiclassical orbits in the plane {fj,,9p): 
fJ-i{Op) = -£/w + [ri,'k±] ■ Zq. 

The quasiclassical orbit in the (/i, 6'p)-plane for a sin- 
gle Abrikosov vortex is shown in Fig. [Ijb). Each point 




FIG. 1: Schematic plot of trajectory precession around vortex 
line in real space (a) and corresponding quasiclassical orbit for 
e = in the (/i, ^p)-plane (b). Vortex core is shown by grey 
circle. 

at this orbit corresponds to a straight trajectory pass- 
ing through the vortex core [Fig. [Ija)]. Precession of 
the quasiclassical trajectory is described by the Hamilton 
equation: hd9p/dt = de/dfj, which provides us the pre- 
cession frequency fl = —uj/h. This precession is a result 



of the small deviation from the exact Andreev backscat- 
tering of quasiparticlcs in the vortex core. In Fig. [Hb) 
the direction of the trajectory precession is shown by 
arrows. The discrete spectrum of subgap quasiparticle 
states can be found using the Bohr-Sommerfeld rule. In 
the case of several vortices we have several crossing qua- 
siclassical orbits in the (/^, 0p)~plane. These orbits are 
shown by dash lines in Fig. HJa) for a particular case 
of two vortices with ri = (—a/2,0) and r2 = (a/2,0). 
Each crossing point of quasiclassical orbits fJ-i{9p) and 




FIG. 2: Schematic plots of quasiclassical orbits for e = in the 
(/i, 6'p)-plane (solid lines) for (a) two vortices with intervortex 
distance a and (b) vortex near the flat surface (a/2 is the 
distance between vortex and surface) . The orbits for two non- 
interacting vortices are shown by dash lines. 

fJ-j{9p) corresponds to the trajectories passing through 
the cores of i—th and j— th vortices. It is natural to ex- 
pect that the degeneracy at these points will be removed 
if we take account of a finite probability of quasiparticle 
tunneling between the cores. Let us consider the vicin- 
ity of the degeneracy point, e.g., 9p ~ [see Fig. [2Ja)]. 
The trajectory characterized by the angle |6'p| ^ £,/a 
passes through both vortex cores, and therefore the wave 
function along such trajectories can be found as a su- 
perposition of two states localized at different vortices 
and having close energies: Syi = — 'j-'[m ^ (fcj^a/2) sin^p] 
and £1,2 = —uj[fi+ {k±a/2) sin^p]. The transformation of 
the quasiclassical spectrum occurs due to overlapping of 
the corresponding wave functions can be described using 
a standard quantum mechanical approach describing a 
two-level systemic, which yields the secular equation 



(e-e„i)(e-e„2) = {Sef 



(3) 



and resulting splitting of isoenergetic lines near the de- 
generacy point {9p ~ for our example): 



-ujfi±^Lu^kA^a/2)^9l + {Se 



(4) 



The tunneling of quasiparticlcs between vortex cores is 
determined by the exponentially small overlapping of 
wave functions localized near the cores and results in the 
splitting of energy levels: 6e ^ Aoexp{—kFaij/{k±^)), 
where a^- = [r^ — r^l is the distance between vortex 
lines, and kpaij/k± is the distance between vortex cen- 
ters along the trajectory. The estimate for the splitting 



(5/i ~ fc/w of isoencrgctic lines in the (^, 9p)-plane reads 

[seeEq. (HI)]: 



Sfi{aij) ^ fcj_^cxp 






(5) 



As a result, we get the orbits fj,*{9p) with a qualitatively 
new behavior [solid lines in Fig. [H^a)]: each of these or- 
bits consists of parts corresponding to the classical quasi- 
particle trajectories passing through the cores of different 
vortices. 

The tunneling between the cores of the different vor- 
tices becomes significant when the energy splitting Se 
is comparable with the interlevel spacing wq, i.e., when 
(5/U > 1 in Eq. (O. According to the above condition on 
Sfj,{aij) the tunneling is most efhcient for k± = kp and 
o-ij < cic, where Oc — £,hi(kFO i^ ^ critical intcrvortex 
distance. Using the percolation theory language we can 
consider the vortices to be bonded if a-y < Oc and de- 
fine a cluster in a flux line system as a set of M vortices 
bonded either directly or via other vortices. Certainly in 
mesoscopic superconductors the cluster dimensions i.„ 
can not exceed the sample size. The cluster is charac- 
terized by a set of hybridized quasiparticle states: with a 
change in the kj^ -direction the wave function experiences 
a number of subsequent transitions between the cores of 
neighboring vortices. Taking, e.g., the upper quasiclassi- 
cal orbit in Fig. [HJa) we obtain the wave function con- 
centrated near the cores of the right and left vortices for 
the angular intervals < 9p < n and tt < 0p < 27r, re- 
spectively. Further decrease in the intervortex distance 
results in the increase in the tunneling probability and, 
thus, the increase in Sfi{aij). Finally, for a^ —^ we get 
a set of M lines fi = const parallel to the 9p axis, i.e., M 
anomalous branches crossing zero energy at angular in- 
dependent impact parameters and corresponding to the 
Af -quantum vortex. Certainly this limit can be realized 
only in mesoscopic samples. 

Within the quasiclassical approach one can estimate 
the intervortex tunneling efficiency using the Landau- 
Zehncr transition theory. Let us consider the vicinity of 
the degeneracy point, e.g., 9p ^ [see Fig. [D^a)]. The 
tunneling probability of transition from one quasiclassical 
orbit to another is given by the cxpressioni^ 



M^ = exp I -4Ini / n{9p) d9p 



(6) 



where 9* = 2de/{ujk±a) and fJ.{9p) should be taken from 
the Eq. ^ with lower sign. Finally, we obtain the fol- 
lowing estimate for the tunneling probability: 



W = exp(-27r(5^/A^)2) 



(7) 



where A/i = \/k±a is the quantum mechanical uncer- 
tainty of the angular momentum. Therefore, we can ne- 
glect the tunneling between quasiclassical orbits while 
S^i > A/i. 



Following Ref. [l7| one can obtain the discrete energy 
levels applying the Bohr-Sommerfeld quantization rule 
for canonically conjugate variables fi and 9p: 



ld{9p)d9p = 2Tr{n + (3), 



(8) 



where n and ng are integers, 2Tmg is the period of the 
fJ-{9p) function {1 < ng < M), and /? is of the order unity. 
The period of /i(0p) can be larger than 27r {rig > 1) if the 
Landau-Zehner transitions between quasiclassical orbits 
are not negligible. Depending on the ratio Sfi/A/j, one 
should apply this quantization rule either to the orbits 
fJ'i{9p) or to the orbits ^*{9p). In the momentum region 



1 - [min(a.y)/ac]2 < |fe| < kp 



(9) 



we can neglect the splitting of isoenergetic lines [5^ <C 
A/i] and Eq. ([5]) written for the orbits i-ii{9p) gives us the 
CdGM spectrum with a minigap ujo/2 ~ cj(fc^ = 0)/2. 
For min(aij) > Qc the CdGM expression holds for the 
entire momentum range. For vortices forming a cluster 
the quasiparticle states bonded by intervortex tunneling 
appear in a finite momentum interval \kz\ < k*, where 






-yl - [min(ay )/ac]^ 



(10) 



In this limit the quasiparticle tunneling between the cores 
results in the qualitative modification of spectrum which 
can be obtained by substituting fJ.*{9p) into Eq. dS]): 



i("^2 



Ao 



-/3 



+ h{ri,..rM) 



(11) 



where i — I..AI. The spectrum (jlip is similar to the one 
of a multiquanta vovteii^'^^ which recovers in the limit 
fly —^ when \bi\ < ^. The multi-vortex cluster geome- 
try and its size L^, determine the effective impact param- 
eters bi{ri, ..rj\/) which vary in the range — i„ < 6,; < Ly. 
Taking a two- (three-) vortex cluster with ^ < a < Uc 
as an example we get 61,2 ~ ±a (61,3 ~ ±a, 62 = 0). 
Contrary to the CdGM case the spectrum branches (fTTI) 
can cross the Fermi level as functions of k^ as we de- 
crease the intervortex distance a and the minigap is sup- 
pressed. The DOS consists of M sets of van Hove singu- 
larities corresponding to the extrema of Eniikz) branches. 
The energy interval between the peaks belonging to each 
set is ojo- For a fixed energy the DOS as a function of 
a exhibits oscillations with the period of the order of 
the atomic length scale. Experimentally the intervortex 
distance can be controll ed by a varying magnetic field. 
For typical values a ~ y^(j)o/H we get the following field 

scale of DOS oscillations: SH/H ~ y^hujH/sp, where 
£f is the Fermi energy. The oscillatory behavior should 
affect both thermodynamic and transport properties at 
low temperatures though in real experimental conditions 
the DOS peak structure is certainly smeared due to the 



various mechanisms of level broadening, e.g., finite tem- 
perature, fluctuations in vortex positions, impurity scat- 
tering effects, etc. It should be noted that for typical 
values kp£^ = 10^ — 10'^ the critical distance Uc/^ ~ 4 — 6 
exceeds the core radius and the spectrum transformation 
starts at the fields H ^ 4'Q/a'^ ^ ifc2[ln(fcF0]~^ when 
the vortices are indeed well-separated. 



B. Effect of normal reflection of quasiparticles at 
the boundary 

The above discussion of the behavior of anomalous 
branches assumed a negligible role of the normal scat- 
tering of quasiparticles at the sample boundary. Such 
assumption is certainly no more valid when the intervor- 
tex distance becomes so large that some of the individual 
vortices approach the sample boundary and their spec- 
trum is determined by the interplay of Andrcev reflection 
and normal scattering at the boundary. For rather small 
samples this interplay can influence the spectrum even 
for a vortex positioned at the sample center— li^. 

In order to focus on the role of the boundary effects we 
consider a model situation when the intervortex quasi- 
particle tunneling can be neglected (i.e. a.y > Oc) be- 
cause of the rather large cluster size. Thus, we can study 
the spectrum modification for a single vortex approach- 
ing the boundary characterized by a certain curvature 
in the plane perpendicular to the vortex axis (see Fig. 
131) . For the sake of simplicity we restrict ourselves to the 
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FIG. 3: The geometrical optics analogy for the reflection of 
quasiclassical trajectories at the concave {F < 0) parabolic 
boundary (left panel) and the corresponding modification of 
the isoenergetic lines (right panel) for d > —F (a), (c), and 
for d<-F (b), (d). 

case of a smooth and specularly reflecting surface. Ob- 



viously the quasiclassical spectrum should be disturbed 
most strongly for trajectories which return back to the 
vortex core after the normal reflection at the boundary. 
These trajectories experience the reflection from a rather 
narrow surface region near the point positioned at the 
surface at a minimal distance d from the vortex center. 
In this region we may consider the surface profile as a 
parabolic cylinder with a certain focal distance F. In- 
troducing a polar coordinate system (r, 0) with the ori- 
gin in the vortex center one obtains the following equa- 
tion describing this parabolic cylinder at small 6 angles: 
r{e) = d{l + [1/2 + d/{'iF)]e^}. Note that we should put 
d/F > —2 so that the distance d is indeed a minimal dis- 
tance to the surface. As it is shown in Fig. [3] the vortex 
center is positioned at the optical axis of the parabolic 
mirror. The trajectories experiencing normal reflection 
and passing twice through the core can be considered 
within the paraxial approximation. 

In this case the scattering rules for trajectories reflect- 
ing from the surface can be obtained employing an anal- 
ogy with a textbook picture describing the system of rays 
and images in geometrical optics. For a particular case 
of a concave parabolic mirror {F < 0) the system of rays 
is schematically shown in Figs. EJa), (b). For the ob- 
ject (white solid arrow) situated at the distance d from 
the mirror the image (white dash arrow) is formed by 
the reflected rays at the coordinate / = —d/h, where 
h = — (1 + d/F). The type of image is determined by the 
sign of /i: if /i > the image is real, i.e. it has the coordi- 
nate / < and is situated at the same side of the mirror 
relative to the object [see Fig. [SJa)], otherwise the image 
is virtual with / > and situated at the other side of the 
mirror [see Fig. [3[b)]. Using these well-known results 
the parameters of reflected trajectories can be derived 
from the simple trigonometry. 

Let us consider the trajectory which makes a small 
angle |6'p| <C 1 with the a;-axis and has a small impact 
parameter \b\ ^ d. Then, the reflected trajectory has the 
angle 9p — ii + hOp and its impact parameter is 6 = hh. 
The impact parameters of incident and reflected trajec- 
tories are deflned relative to the point = Q, r = d 
positioned at the surface at the minimum distance to the 
vortex center. This point coincides with the coordinate 
system origin in Figs. [31Ia),(b). 

In the (/i, 0p)-plane one can define isoenergetic lines 
corresponding to the incident and reflected trajectories. 
The intersection of these lines occurring for \0p\ <?C 1 cor- 
responds to the situation when both the incident and 
reflected trajectories pass through the vortex core. The 
degeneracy at the intersection point should be removed 
by the splitting of isoenergetic lines due to the interaction 
of vortex core states. One can obtain two qualitatively 
different regimes of splitting, determined by the ratio be- 
tween the focal distance F of the parabolic mirror and 
the distance from the vortex to the surface d, i.e. by 
the sign of h. The splitting corresponds to the continu- 
ous transition from one isoenergetic line to another. The 
velocity of motion along isoenergetic lines is determined 



by the angular velocities of trajectory precession, given 
by ri = dOp/dt and Cl = dOp/dt = Ml for the incident 
and reflected trajectories, correspondingly. Therefore, if 
Q, and f2 have the same signs (/i > 0), then the directions 
of precession of the incident and reflected trajectories co- 
incide. In this case, the orbit transformation is analogous 
to the one. that have been obtained for the pair of inter- 
acting vortices [Figs. [3Jc) and[2)Ja)]. Otherwise, if /i < 0, 
the incident and reflected trajectories precess in differ- 
ent directions, therefore the splitting occurs as shown in 
Figs. WiA) and[2{b) analogously to the transformation of 
isoenergetic lines which can be obtained for the interact- 
ing vortex and antivortex. 

To study the spectrum transformation let us consider 
in detail the simplest configuration: the vortex line is 
situated at the point (—a/2, 0) near the flat boundary 
of superconductor, occupying the half-space a; < 0. If 
we neglect the normal reflection of quasiparticles at the 
boundary, the isoenergetic line in the (/i, 6'p)-space is 
given by iiv(9p) = —e/uj — {k±a/2)sm.9p, shown for 
the particular case e = by the dash line in Fig. 
[5Jb). The flat boundary is characterized by the infi- 
nite focal distance: F = oo, therefore h = —1 and 
we obtain the mapping of incident and reflected tra- 
jectories according to the following rule: b — —b and 
Op = TT — 9p. Then, we obtain another isoenergetic line: 
l^aviPp) ~ — yUtiCTT ~ Qp) = e/w -I- (fcj_a/2) sinSp, which 
corresponds to the reflected trajectories and is shown in 
Fig. lUb) by another dash curve. Note that the isoener- 
getic line fiavi&p) (with the opposite direction of trajec- 
tory precession) coincides with the isoenergetic line cor- 
responding to an antivortex placed at the point (0, a/2) 
outside the superconductor. It means that the spec- 
trum of the vortex near the flat surface can be obtained 
considering the spectrum of the vortex-antivortex sys- 
tem. Indeed, the pair potential of the vortex-antivortex 
system is invariant under the reflection at the x ~ 
plane: A(x,y) = A{—x,y), yielding the symmetry of 
the quasiparticle wave function: ^(x,y) = ±^(— x,j/). 
The odd wave functions obey the boundary condition 
^{Q,y) = 0. The even wave functions obey the bound- 
ary conditions d'i!{0,y)/dx = and the corresponding 
energy levels should be omitted in order to obtain the 
spectrum of the vortex near the surface. The isoener- 
getic lines ^v.av{Op) = ±[e/a; + (fc_La/2) sin^p] intersect 
at certain points, e.g., at Op = im for e = 0, where n is in- 
teger. The degeneracy at these points is removed by the 
splitting of isoenergetic lines, shown in Fig. [2jb) by the 
solid lines. Considering the interaction of the two quan- 
tum states with close energies e„ = —uj[iJ,+ {k±a/2) sin 6'p] 
and Eav = LL>[^—{k±a/2) sin 6'p], we obtain again the secu- 
lar equation ([3]) with Sui = £v and £„2 = £av Therefore, 
the quasiclassical orbits near the degeneracy point Op ~ 



are given in this case by the following expression: 

e = ±^iiJ^iy + iSe)^- uj{k^a/2)9p . 



(12) 



The classically forbidden angular domain at e = has 
the width SOp — 4:Se/{u!k±a). One can assume that the 



appearance of such classically forbidden domain explains 
the deep structure in the local DOS profile observed nu- 
merically in Ref. [l3 for a vortex near the flat boundary of 
an s-wave superconductor. As we show below, the classi- 
cally forbidden angular domain results in the suppression 
of the overall DOS and we propose that this mechanism 
should be responsible for the anomalous spectrum branch 
disappearance when the vortex exits the sample. 

To estimate the tunneling probability between the qua- 
siclassical orbits we again apply the theory of Landau- 
Zehner transitions. The Landau-Zchner tunneling prob- 
ability is expressed as follows: 



W = exp -4Im 



9p{fi)dfi 



(13) 



where fi* — Se/uj and Op{^) should be taken from Eq. 
P^ with the upper sign. Finally, we obtain the tun- 
nehng probability as W ^ exp{—2Tr{69p/ AOp)"^) , where 
AOp ~ {k±a)~^''^ is the quantum mechanical uncertainty 
of the trajectory orientation angle. Thus, we can neglect 
Landau-Zehner effects while SOp > AOp, i.e., for a < ac, 
where ac ~ ^ln(A;_L^) is the critical distance, which ap- 
pears to be the same as for a two-vortex system. 



III. ANDREEV EQUATIONS 

Our further consideration is based on the Bogoliubov- 
de Gennes (BdG) equations for particle- (u) and hole- 
like (v) parts of the wave function, which have the fol- 
lowing form: 

- 7— (V^ 4- kl.) u + A{v)v = su , 
Zm 

-— N"^ + kl) V + A* {v)u = ev. (14) 

2m 

Here A(r) is the gap function and r = {x,y) is a radius 
vector in the plane perpendicular to the magnetic field di- 
rection. We assume the system to be homogeneous along 
the z~ axis, thus, the k^- projection of the momentum is 
conserved and the wave function takes the form: 



{u,v) 



ik-,-^ 



^(r) 



Then, the two-component wave function ^ in the mo- 
mentum representation can be written as follows: 



*« = (2iF//I^"^''^^''p- 



(15) 



Let us introduce the polar coordinate system in momen- 
tum space p = p (cos^p, sin^p) = ppo- Then, the coor- 
dinate operator can be written as follows: 



., d ., / d i , , , 

r = ih—- = ih poT^ + - Zo,Po yU 
ap \ op p 



where operator fl of z-projection of angular momentum 
is given by the expression: 



" ![ 1 ■ ^ 



(16) 



Next, we assume that the quasiparticle wave function 
can be viewed as a wave packet with momenta abso- 
lute values close to hk±. This assumption is valid with 
very good accuracy in most superconductors, since the 
characteristic length scale of envelopes of quasiparticle 
waves is determined by the superconducting coherence 
length ^, which is typically much larger than the quasi- 
particle wave length {kp£^ ^ 1). Therefore, one can put 
p = hk± + q {\q\ <C hk±) and obtain: 



ihpo 



d_ 
dq 



I 
2k\ 



[zo,PoJ 



_5_ 
9^ 



where {...} is an anticommutator. Of course, such ap- 
proximation is broken for a very small portion of quasi- 
particles which propagate very close to the vortex axis 
[27r/(A;_L^) > 1]. Let us now introduce a Fourier transfor- 
mation: 



1 



+ 00 



V'P = 7- / 'A(s 



*«"/'' ds . 



(17) 



The variable s is a coordinate along a quasiclassical tra- 
jectory, which is a straight line along the direction of the 
quasiparticle momentum. The trajectory orientation an- 
gle is given by the 9p value. The wave function in the 
real space can be found from Eqs. (|15ll7p : 



27r 







^{r,0)^ I e'k^'^^°'^'-'-H{rcos{e-ep),ep)-^, (18) 






where (r, 6, z) is a cylindrical coordinate system. The ex- 
pression for coordinate operator in (s, 0p)-representation 
reads: 

, d 



Then, BdG equations (fT4|) in (s, 6'j,)-representation 
take the form HiIj{s, Op) = eip{s, Op) with the Hamilto- 
nian given by 



H 



-la^- 



h^k^ d 
ds 



A(f) 
A*(f) 



(19) 



where ax, (Jy, ct^ are the Pauli matrices. 

Note that the gap function operator A(f) in Eq. 
(fT9|) contains a differential operator d/dOp, therefore the 
above quasiclassical equations are still rather complicated 
partial differential equations. A further simplification 
can be obtained considering eikonal approximation for 
the angular dependence of wave function: 

^{s,ep) = e'^^^''^^g{s,Op), 



where 

-klWp^'^'^^ 

is an impact parameter of a quasiclassical trajectory. As- 
suming that the angular dependence of g{s, Op) is rather 
slow, one can neglect its angular derivatives in Eq. (J19p . 
The resulting Andrecv equations characterizing the be- 
havior of the wave function along a trajectory with a 
certain orientational angle Op and an impact parameter 
b read: 



/i^fcj_ dg 
m ds 



(7j;ReA(a;, y)g - (TyImA(a;, y)g = eg 



where 



X — s cos Op — b sin Op . 
y = s sin Op -\-b cos Op . 



(20) 



(21) 



Note that the Andreev equations (|20p can be obtained 
directly from the initial BdG equations p^ if one applies 
the coordinate system transformation (PT|) and neglects 
the second-order derivatives of the wave function. 



IV. QUASIPARTICLE SPECTRUM OF A 
MULTIQUANTUM VORTEX 

We start our quantitative analysis of quasiparticle 
spectra with the case of a multiquantum vortex with vor- 
ticity M. In this section we neglect the effect of normal 
quasiparticle scattering at the sample boundary and fo- 
cus on the peculiarities of the spectrum depending on the 
vorticity value. We take the gap profile in the form: 



A = Dm (r) e^''' 
In s, Op variables one obtains: 



A = D 



M [Vs^ + b^) e' 



Me„ 



s + ib 
Vs2 + 62 



M 



(22) 



(23) 



Due to the cylindrical symmetry the Op- dependence of 
the function g can be excluded using the gauge transfor- 
mation 



g = exp{iMa,0p/2)f 
and quasiclassical equations (|20p take the form 
n^k_L df 



la. 



m ds 
Here we introduce the functions 



+ ^xGr/ - ayGif = ef 



(24) 



(25) 



Gr = Dm 
Gi = Dm 



52) 



Re 



g2 j^ ^2 j Jj^ 



s + ib 


[Vs^ + b^\ 


s + ib 


IVs^ + b^l 



M 



M 



To apply the method analogous to the one used in Ref. 
|8| for a singly-quantized vortex we note that the exact 
solutions of the above equations corresponding to £ = 
can be found in case G/ = 0: 



/± = (l,±i)cxp ± 



'h?kj_ 



Gads 



(26) 



Provided G/? is an odd function of s, which tends to a 
certain nonzero value for |s| — > cc one of these solutions 
appears to decay both at negative and positive s and, 
thus, we get a localized wave function corresponding to 
a midgap bound state. Using this localized solution as a 
zero-order approximation for the wave function the spec- 
trum can be found within the first order perturbation 
theory assuming that |£| ^ Aq. 

For an arbitrary value of vorticity the function Gr 
is not necessary odd. In order to use the perturbation 
method described above we apply a gauge transforma- 
tion 



/ = 



i(7.& 



Vs2 + 62 



(27) 



so that new wave function w satisfies the following equa- 
tion: 



m ds 



^xGj^ 



^yGf' 



£d 



w 



ew 



Here 



£d = -- 



fi^k^ ah 



52 _|_52 



(28) 



(29) 



is the Doppler shift, and 
G^"' = Dm (^s^+^) Re ■ 



Gf ) = Dm ( Vs2 



62^ 



Im ■ 



s -\- ib 


Vs2 + 62j 


s + ib 


^/s^ + bA 



M-2a 



M-2o 



,(30) 



are the real and imaginary part of the off-diagonal poten- 
tial, correspondingly. One can see that choosing M — 2a 
to be an odd positive integer we can change the parity of 
the potentials in the Hamiltonian so that G^ (s) is an 
odd function. Comparing Eq. (^5]) and Eq. (^5]) we ob- 
serve that the above gauge transformation also produces 
a Doppler shift of the energy levels. In principle, both 
the Doppler shift term and the term o-yGj are not small 
and can be of the order of Aq if the impact parameter is 
rather large: 6 ~ ^. Thus, strictly speaking we can con- 
sider the expression ([26l) with Gr replaced by G}^ (s) 
as a zero-order approximation and use the perturbation 
technique discussed above only provided the energy cor- 
rections arising from the terms ayGj and Sd almost 
compensate each other. The anomalous branches which 
we study cross the Fermi level at certain impact param- 
eters —fij/kj_, thus, the perturbation method should be 



adequate in the vicinity of these points. Changing a in 
the interval < a < M/2 we get a set of possible odd 
M — 2a values providing us a set of different zero-order 
approximations which allow to obtain the spectrum as a 
function of b. Surprisingly, we shall see below that this 
method can describe the spectrum behavior even beyond 
its validity domain, i.e. when the energy is comparable 
to ±Ao. 

It is convenient to parameterize the wave functions as 



=C(s) 



gJr,(s)/2 
„-J77(s)/2 



(31) 



As a result, we have the equations 
ft^fcj^ drj 



2m ds 

h^k± dC 

m ds 



]j. ' cos rj + Gj smi] = e — Ed , 



.(") 



(a) 



G^"^sin77-G}"'cos?7 = 



(32) 



Considering the localized states one should take the fol- 
lowing boundary conditions for odd M — 2a values: 



cos77(±(X)) = ±e/Ao , sin77(±oo) = JI—e^/Aq. 

(33) 
In order to construct the solution we note that the mutual 
phase 77 of the electron and hole components in the zero- 
order solution (|^5|) is constant [77(5) = tt/2], therefore, 
within the perturbation theory we can linearize the Eq. 
([32| for 77 close to tt/2 introducing 7'7 — tt/2 — -q: 



dfj 
'd's 



2m 
Wk 



-G 



R V = 



2m 
Wk7 



£d- 



g\"^ - e 



(34) 



To exclude the divergent solutions of this equation 
we should impose the integral condition describing the 
anomalous spectral branches: 



'-M ~ 



/„°° Sd + G^'is) 



(")/ 



-^(^) ds 



J^ e-^(^) ds 



where 



K{s) 



2m 



h?k 



J- Jo 



G^^\t)dt. 



(35) 



(36) 



Taking Di{r) = Aqt/^t^ -|- ^2 and a = for the sim- 
plest case of a singly-quantized vortex (A/ = 1) we get 
an explicit expression for the CdGM spectrum: 



-(0) 



A„5 /Co (2777AoV^^+eJ/(^'fc±) 
V^' + il /Ci (2777Aov/^^TeJ/(^'fci) 



(37) 



where /C„ is the McDonald function. To obtain the spec- 
trum for M = 2 we need to consider the case a = 1/2 
so that to reduce the problem to the one for a unity vor- 
ticity and a certain Doppler shift. For the case M = 3 
the set of anomalous branches can be obtained taking 



two a values: a = and a — 1. The solution £3 (&) 
gives us the branch crossing the Fermi level at zero im- 
pact parameter, while two other branches are described 
by £3 ' (6) . The typical plots of quasiparticle spectra for 
vortices with winding numbers M = 1,2,3 are shown 
in Fig. m Comparing the spectra e^ with the bran- 
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FIG. 4; (Color online). The anomalous spectral branches as 
functions of the impact parameter h for fcz = 0, obtained 
from Eq. (I35j (blue dash lines) for M = 1 (a), M = 2 
(b), Af = 3 (c). The anomalous spectral branches obtained 
from numerical solution of the eigen-value problem (|14|l are 
shown by solid red lines. T he gap profiles are approximated as 
DA/(r) = A()(r-/i/r2+^2)*^ with the parameter kpi = 200. 



ches obtained from our direct numerical analysis of BdG 
equations ()14p one can see that our perturbation method 
provides a reasonable description of the low energy spec- 
trum behavior. The numerical solution of the eigen-value 
problem P^ was carried out using a representation of the 
BdG operator in the truncated basis of the normal-metal 
eigenfunctions. We solve the system of linear equations 
which correspond to the boundary conditions at the su- 
perconductor/insulator boundary of a cylinder with a fi- 
nite radius (we took R — 7^). The small oscillations of 
the spectrum as a function of h (solid red lines in the Fig. 
|4]) result from the interference of incident and reflected 
from the boundary quasiparticle waves^^. Note, that our 
perturbation procedure fails to describe proper behavior 
of all the branches in the vicinity of the gap value Aq: 
e.g., for M — 2 the function e\ jumps at fe = from 
the upper branch to the lower one and, thus, we can not 
describe the part of the upper branch approaching Aq for 



6> 0. 

To obtain the spectrum as a function of discrete vari- 
able /i instead of a continuous h we should apply the 
Bohr-Sommcrfcld quantization rule for the angular mo- 
mentum [see Eq. (|5])] with /? = {M/2} which results from 
the obvious condition that the wave function g is single- 
valued. Here {...} denotes the fractional part. For the 
odd (even) vorticity M we obtain /j, = n -I- 1/2 {^ = n), 
where n is an integer. 



V. DECAY OF A MULTIQUANTUM VORTEX 
INTO A SET OF SEPARATED VORTICES 

In this section we consider modification of the quasi- 
particle spectrum caused by the decay of a multiquantum 
vortex into a set of separated singly-quantized vortices, 
which occurs under the magnetic field decreasing. For 
simplicity sake wc restrict ourselves to the case of two- 
vortex system with a certain intervortcx distance a con- 
trolled by the external magnetic field. The case a — 
corresponds to a doubly-quantized vortex while the limit 
a ^ ^ corresponds to a pair of isolated singly-quantized 
vortices. In this section we again neglect the effect of 
normal quasiparticle scattering at the sample boundary. 



A. Quasiclassical consideration 

As a first step of our analysis we choose to apply the 
approximate quasiclassical procedure developed in the 
previous section. It is natural to expect that the va- 
lidity range for this method should be restricted to the 
region of rather small distances a < Gc when one can 
neglect the Landau-Zener tunneling between the quasi- 
classical orbits fJ.{9p) described in Section[TIl To describe 
the system of two singly-quantized vortices positioned at 
r = ±a/2 = ±(a/2,0) wc fit the gap function as follows: 



A(r) = Ao/i(r-a/2)/i(r + a/2) 
= Ao|/i(r-a/2)||/i(r + a/2)| 
X + iy — a/2 x + iy + a/2 
\x + iy — a/2\ \x + iy + a/2\ 



(38) 



where /i(r) is a normalized gap function of a singly- 
quantized vortex. It is convenient to rewrite the above 
expression as a superposition of functions with two dif- 
ferent vorticities: 



A(r) = Ao[/2(r)e2'^ + /o(r)] 
Taking the simplest core model 



(39) 



(40) 
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x^ + y"^ 



with the core size ^„ we obtain: 
/2(r) = 
/o(r) = - 



aV4 

v/(a;2 + 2/2 + ^2 _|_ 0,2/4)2 _ ^23,2 ■ 



To solve the equations ((20|) with the gap function given 
by p9p we apply the gauge transformation (|27p with a = 
1/2. The expression for the quasiclassical spectrum takes 
the form (1351) with 



Gr = An 



h{x,y) 



h{x,v) 



Vs2 + 62 

scos(26'p) — 6sin(26'p) 



G/ = Ao 



h{x,y) 



Vs2 + 62 
6 



fo{x,y) 



Vs2 + 62 
ssin(26'p) + 6cos(26'p) 
Vs2 + 62 



(41) 



(42) 



where (x,?/) variables are given by Eqs. (PT|) . 

The resulting dependencies of the impact parameter vs 
9p for zero energy and different intervortcx distances are 
shown in Fig. [5l These calculations of the low energy 
part of the spectrum appear to be in a good agreement 
with the two-level model ^ and, thus, the angular de- 
pendence of the impact parameter can be fitted by the 
expression; 



HOp 



uj{kz,a)k^ 



±\bHh,< 



a . 
— sin I 
2 



(43) 



where (D(fe2,a) ~ Ao/fc_L^, and 6 = d£/{k±uj) is the 
splitting of quasiclassical orbits. In the limit a = 
we get the spectrum of a doubly-quantized vortex: s = 
oj{kz, 0)[6±6(fcz, 0)], where 6(fc^, 0) is of the order of ^ for 
small kz values. For large intervortcx distances the value 
6 is exponentially small: 6 ~ ^ ex.p{—kFa/ {k±^)) [see Eq. 
(ini)]. Generally, in the whole interval of distances a the 
splitting of quasiclassical orbits is defined by the over- 
lapping of the wave functions localized in the cores of 
neighboring singly-quantized vortices. This overlapping 
is determined by the factor exp{—Ko{s)) describing the 
decay of the wave function in a singly-quantized vortex 
(see Ref. i): 



Koi-s) 



2m 



h^k 



-L Jo 



A(i) dt . 



(44) 



Thus, the spectrum p5|) can be fitted by the 
one describing a two-level system if we put 6 = 
6(fcz,0) exp(— iiro(a/2)). Taking the vortex core model 
(itol) we obtain: 



6 = 6(fc^,0)exp -2 



kp £.v 




(45) 
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FIG. 5: The isoenergetic curves b(Sp) for a two-vortex system. 
We choose here e = 0, fcz = 0, a = 0.3^ (a), a = C (b), a = 3^ 
(c). The dash curves correspond to the isoenergetic lines for 
two non-interacting singly-quantized vortices. The gap is 
approximated by Eqs. p8)l . (|40)l with 5„ = ^. 



Comparing this approximate expression with the orbit 
splitting calculated using the spectrum (|35p (see Fig. [5] 
for the particular case k^ = 0) we can find appropriate 
6(fc2,0) values. 

To find the quantized energy levels we can use the 
Bohr-Sommerfeld quantization rule (l8|) which takes the 
form S{e,kz) — 27r(n -f /?), where 



S{e,kz) 



2tt 



fi{9p) dOp 



2-K 



b{0p) ddp 



(46) 



is the area under the isoenergetic line /^(er, k^) in the 
(/i, 0p)-plane. Calculating this integral for the two-level 
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FIG. 6: Splitting of quasiclassical orbits b obtained from Eq. 
(|45|l (solid line) and found from the spectrum of a two-vortex 
system pSfl (circles) for kz — 0. The gap profile is approxi- 
mated by Eq. (|39l) with C« = ?■ 



model ([13]) we obtain: 



S{s, K) = -27r- ± 2fc_L V a^ + 462 E ^ , 

(47) 
where E(fc) — J^ v 1 — fc^ sin^ 9 dO is the complete el- 
liptic integral of the second type. As a result, we find the 
quasiparticle spectrum of a two-vortex system 



-n-P± 



k^^/a^ + 462 



E 



Va2 +452 



(48) 

The spectrum (|48|) is analogous to the spectrum of the 
doubly-quantized vortex ([1]) and is drastically different 
from the CdGM spectrum. The case a = gives us the 
spectrum ([T]) of a doubly-quantized vortex. Treating the 
opposite limit a 3> 6 with the logarithmic accuracy we 
obtain: 



S{e,K 



and 



-27r- ±2fc_La 



In 



(49) 



-71 ± 



k±a 




(50) 



B. Landau-Zehner tunneling between quasiclassical 
orbits. 



The above description is valid provided the intervor- 
tex distance is rather small (a < Qc) when the splitting 
5^ between isoenergetic lines ^J,{Op) is large compared to 
the quantum mechanical uncertainty of the angular mo- 
mentum A^ and the probability of Landau-Zehner tun- 
neling between quasiclassical orbits can be neglected [see 



Eq.([7])]. As a next step, we proceed with a quantitative 
analysis of quasiparticle spectrum of the two-vortex sys- 
tem in case of the small splitting of the isoenergetic lines 
in (/i, 0p)-plane, when the vortices are well separated 
so that a ^ ac- In order to calculate the quasiparticle 
spectrum taking into account the influence of Landau- 
Zehner tunneling we should go beyond the quasiclassical 
consideration of the angular precession of quasiparticle 
trajectories. It means that we can not neglect the non- 
commutativity of canonical variables: [fl,Op\ = —i, where 
9p is the trajectory orientation angle and fi ~ —id/dOp 
is the angular momentum operator. Keeping in mind 
the symmetry of the gap function in a two-vortex sys- 
tem we can reduce the problem to the one describing a 
single vortex with an additional boundary condition im- 
posed on the wave function at the plane x = positioned 
between vortices. Indeed, the gap function distribution 
corresponding to the two-vortex system possesses the fol- 
lowing symmetry: A(x,y) = A*(— a;,j/). As a result, for 
the eigenfunctions we obtain: 



^{x,y) = e'^i!*{~x,y), 



(51) 



where x is a constant phase. The spectrum does not 
depend on x, since for any eigenfunction ^ satisfying 
Eq. ([ST]) we can introduce a new function ^i = ^e~''<^/2^ 
which corresponds to the same energy level and has the 
following symmetry: ^i(x,?/) = \E'|(— cc,?/). Therefore 
we can choose x = ^^^ obtain the boundary conditions 
at the plane x = Q: 



^ :^'^* 



dx 



dx 



(52) 



For the sake of simplicity we neglect the anisotropy of the 
gap function around the vortex positioned in the half- 
space X < Q. Nevertheless the solution can not be char- 
acterized by a definite angular momentum because of the 
above boundary condition responsible for interaction of 
different angular harmonics. Thus, following Ref. [ij we 
introduce the angular-momentum expansion for the so- 
lution: 



^{s.Op) 






^ifj,Bp+icT^8p/2 



5m(s), 



(53) 



where ^ = n + 1/2, and n is an integer. The function 
(jf^ satisfies the Andrecv equation pO|) along the quasi- 
classical trajectory with 6 = ^fj./k±. For small impact 
parameters 6 <C ^ the Eq. ([^0]) can be solved analytically, 
yielding a general expression for the function 5^(5) in the 
following form: 

where c^, d^ are arbitrary constants and we choose 
the fundamental solutions so that Gi^(— 00) = while 
G2^(+oo) == (see Ref. [H: 



G2fi ~ 



e-l^"(^)l/2-i^(sgms 



DazC 



\Kois)\/2 



-|Ko(s)l/2 _ ,-2 



i-(sgns-l)<7,el^«(")l/2 
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Here A = (cxp(i7r/4),exp(— i7r/4)), 



A 



^0 



k±Uo 



(54) 



(55) 



Here the CdGM spectrum is taken as a linear function 
of \fj,\ for small /j <^ fc±Cti- ^(p-) ~ — /icj with interlevel 
spacing 

A fc2^e 7o s 

It is convenient to introduce the angle-dependent func- 
tions: 



C{Op) 






D{ej,) = Y,e'>^^''+-U^. 



These functions appear to be nonzero only in the angular 
interval — 7r/2 < 9p < tt/2 because of the decay of the 
wave function ^(x,?/) in the left half-space far from the 
vortex. At the boundaries of this angular interval we 
should impose the conditions 



C(±7r/2) = ±£1(^71-/2) 



(56) 



which result from the continuity of the wave function 
fp{s,Op). 

Within the large angle domain 16*^1 ^ ^/a the wave 
function ^/'(s, 9p) can be found using a tight binding ap- 
proximation as a sum of two single-vortex solutions lo- 
calized on vortices at r = ±(a/2,0) (see Ref. 0). Com- 
paring the tight-binding solution with Eq. (|53p we obtain 
that E7(M)e'^''^c^ = and J2l{tJ-)e''^^''^''^''df, = 0, or 



Ciep),D{9p^n)^e 



— ie6p/uj 



(57) 



The deviations of angular functions C{9p), D{9p) from 
([57)1 due to the intervortex quasiparticle tunneling oc- 
cur in the narrow angular domain \9p\ <^ ^/a. Within 
this domain we apply the stationary phase method to 
handle the integral in Eq. pH)) . For a given value of an- 
gular momentum fi the stationary phase points are given 
by sin(6'p ~ 9) = fi/[kjLR{9)], where R{9) = a/{2cos9) 
corresponds to the line a; = in the polar coordinate 
system with the origin in the vortex center. Assuming 
|/i| ^ k±a and 16*1 <C tt, we obtain the stationary phase 
points: 9pi fn 9 + 2fi/{k±a) and 9p2 = tt + 9 — 2fj,/(k±a). 
Then, the expression for the wave function ^(r, 0) at 
r = R{9) and \9\ ^ tt reads as follows: 



where ip = k±a{l + 0^/2)/2 and the discreteness of an- 
gular momentum /.i is neglected. Then, from Eqs. ([F 
we obtain: 



Im 



1^ /h±a -Ko{a/2) 



i-fa^ ) Xe''"' Cfj,d^i 



= 0, 



which yields the following equation (see Appendix E| : 
(e-^«("/2)-»a,7)AC2(0) 



(58) 



(• 



e-2"^(e-^«('^/2)-hia^7) A*Ci(6l), (59) 



where <f = /c_La(l - 6lV2)/2, C2{9) ^ iC{9) and Ci{9) 
C*{—9). Then, Eq. ((59)) can be written as follows: 



lUJ 



ILO 



d_ 
09 
d_ 
09 



where the overlap integral is 



£JC2{9) = e-^'^JCii9), 
e]Ci{9) = e^'^JC2{9), (60) 



(61) 



7 - ^p-Koi^m 



Analogously, for the functions D2{9) = iD{9) and 
Di{9) = D*{-9) we obtain: 

d 



-e\D2{9) = e''^JD,{9), 
-^Di{9) = e-^'^JD2{9). (62) 



961 

d_ 

'd9 



Considering the Eqs. (|60p. ([62]) in the angular domain 
(^/a) ^ \9\ 3> {£,/a)J/Ao we can neglect the rapidly 
oscillating right hand side, therefore the asymptotic form 
of angular functions C and D corresponds to the non- 
interacting vortices: 



C 
D 



(Ci,C2)-e— ''/"(c~i,£2), 



iDi,D2 



-iee/u 



{di,d2 



(63) 



where ci, C2, di, d2 are arbitrary constants. Then, solv- 
ing Eqs. (pO|) . ([m we can find the transfer matrix X 
matching the large-angle asymptotics: 



7t/l, 



= XC{-9th) , 



*(i?(6l), 9) = e'f / e''^'/'^^" L-Ko(a/2) _ -^^ \ Xe't'Ocpdn 



where X is a transfer matrix, and 9th ^ S,/a is a threshold 
angle where we match the solutions for the large-angle 
and small-angle domains. Introducing new functions 

52(6*) = C2{9)e"^-"''^''^ , 51(6*) = Ci(6l)e-^'^-^^^/'^ , 

we obtain the equations: 

i^ + ~9]b^^pB2, 
o9 I 



-oo 



i-.-9]B2=pBi 



(64) 
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where 



P = 



ujy/k±aj2 



(65) 



and 9 = 9^Jk±aj2. This equations coincide with equa- 
tions obtained in Ref. H- The problem described by the 
Eqs. (|64p is equivalent to the one describing the inter- 
band tunnelingi^ or one-dimensional motion of a Dirac 
particle in a uniform electric field and the solution can 
be written in terms of the parabolic cylinder functions 
(see Appendix IBJ) . yielding the transfer matrix: 



X = e-'^P'/^j _j_ j (o-^Rer -f- (Jxlmr) 



(66) 



Here / is the unity matrix, 



'irp^ /A iX 



T = ^2sinh(7rp2/2)e' 

+ arg(r(l-4^^ - 



(67) 



and r is the gamma function. Analogously, for the func- 
tions T>{0) we find: 

T>{eth) = a^Xa^T>{-eth) ■ 

Matching the wave function in different angular domains 
and using the boundary conditions (|56p we obtain the 
spectrum: 

cos(7r£/w) = ±e-'^p'/4^2 sinh(7rp2/2) sinx ■ (68) 

Within the present theory we can not determine the 
threshold angle 9th precisely. However, since the depen- 
dence of X on 9th is logarithmic, the 6't;i -dependent term 
in Eq. (j67p can be considered as an additional constant 
phase of the oscillations of the energy levels. The spec- 
trum of two vortices calculated using the Eq. (|68p for 
small (a > Oc) and large (a < Oc) intervortex distance is 
shown in Fig. [71 One can see that the transformation 
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FIG. 7: The quasiparticle spectrum of two vortices calculated 
using the Eq. ^ ior a ^ 5^ (a) and a = 3.5^ (b). The 
CdGM spectrum is shown by the dash lines. The vortex core 
profile for a single vortex is approximated by Eq. (|40p with 
C. = C, kpS, = 200. 

of the spectrum e{kz) occurs according to the scenario 



suggested above: as we decrease the distance a below Oc 
the crossover to the doubly-quantized vortex spectrum 
starts in the region of small kz values defined by the con- 
dition p > 1. For p <C 1 we get the CdGM spectrum with 
a small oscillatory correction: 

e~oj[n + 1/2 T (-l)"(p/V^) sin(fc_La + 7r/4)] . (69) 

The minigap £min = [1^2p/y^ciJo/2 vanishes for a = Oc. 
For the large values p S> 1 the Eq. (j68|l yields the 
spectrum in the form: 



n ± 



kj_a + p'^ In (9th\/kj_a/p' 



(70) 



Here we have used the asymptotic formula for the argu- 
ment of gamma function at p 3> 1: arg(r(l — ip'^/2)) « 
— (p^/2)[ln(p^/2) — 1]. This result can be formulated as 
the Bohr-Sommerfeld quantization rule S'(e, fc^) = 2Tm 
with 



S{e,kz 



£ 

27r-± 
a; 



2fcj_a-f 2p2ln 



^ • (71) 




This result coincides with the formula ([49]) obtained by 
the evaluation of the Bohr-Sommerfeld integral at qua- 
siclassical orbits p.(9p) such as ones shown in Fig. [HJb). 



VI. BOUNDARY EFFECTS: VORTEX NEAR 
THE SURFACE 

Now we proceed with a quantitative analysis of the 
effect of normal reflection of quasiparticles at the bound- 
ary. We consider a vortex near a smooth surface approx- 
imated by a parabolic cylinder (see Sec. IHBp . Consid- 
ering the quasiclassical trajectories which are rather far 
from being parallel to the system optical axis we find that 
either incident or reflected trajectory appears to pass far 
from the vortex core. The quasiclassical spectrum for 
this case is the same as for a single vortex: 



e(5, 0p, fcr) = —Lokj_{b — dsin( 



(72) 



This spectrum should be strongly disturbed for trajecto- 
ries passing close to the optical axis when both the in- 
cident and reflected trajectories pass through the vortex 
core. In this case the trajectory orientation angles should 



be taken in the domains \9p\ < ^/d or 



< ^/d and 



the impact parameters defined relative to the point with 
coordinates r ~ d and 9 = (see Sec lHBI for notations) 
should be rather small: |6| ^ d. Solving Eq. (|20p along 
the incident and reflected trajectories and matching the 
solutions to meet the zero boundary condition for the 
wave function at the surface we obtain the quasiclassical 
spectrum for \b\ <C d and \9p\ <C 1 



£{b,9p,kz) 



Ljkj_[il + h)b + 29pd] 



± 



(wfc_L)2(l-/l)252 



J2 . (73) 



14 



This expression for the quasiclassical spectrum can also 
be appUed for the angular interval 1 0p — tt | <^ 1 provided 
we replace the angle Op ^ tt — Op. 
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FIG. 8; The schematic plot of isoenergetic lines fi{8p) corre- 
sponding to the energy e = for a vortex near the parabolic 
boundary positioned at d > |F|: F < 0, h < (a), 
F > 0,h > (b). The solid lines correspond to the trajecto- 
ries passing close to the vortex core. The dash lines represent 
the mapping of solid lines due to the normal reflection of tra- 
jectories at the boundary. 

The Eqs. ([7^ . ([75)1 allow us to determine the form of 
isoenergetic lines n{Op) ~ —k±b[Op). Generally, one can 
distinguish two types of the isoenergetic lines behavior 
near the points Op ~ nn according to the sign of the 
parameter h = -(1 + d/F). If /i < (i.e., if F > or 
F < and d < \F\) there appears a prohibited angular 
domain with the width SOp = p/ \/\p\ [Fig. [H{a)], where 



P 




(74) 



On the other hand, for /i > (i^ < and d > \F\) 
there appears a gap between isoenergetic lines and the 
prohibited angular domain disappears [Fig. Hfb)]. 

Solid curves in Fig. [8] correspond to the trajectories, 
which always pass through the vortex core. The reflection 
of trajectories at the boundary determines the mapping 
of a solid curve on the dash curve. To determine the 
spectrum we should apply the Bohr-Sommcrfeld quan- 
tization rule to the closed path in (^, 6'p)-space corre- 
sponding to the precession of the trajectory around the 
vortex core. Let us study the formation of such closed 
path. We start from the consideration of the trajectory 
on the solid line, precessing along the orbit in the direc- 
tion showed by the arrow in Fig. [51 Starting from the 
point D the trajectory precesses along the solid curve to 
the point A, while the trajectory on the dash curve moves 
from the point B towards the point C. As the incident 
trajectory goes from the point A on the solid curve to the 
point B on the dash curve, the corresponding reflected 
trajectory goes from the point C on the dash curve to 
the point D on the solid curve. Thus, the trajectory is 
reflected from the boundary and appears at the point D 
again, and therefore the quasiclassical orbit in Fig. [5] 
forms a closed path, since the quasiparticlc states in the 



points A and D are identical. Now we can apply the 
Bohr-Sommerfeld quantization rule for the closed orbits 
in Fig. [51 S{e, kz) = 2?™, where n is integer and S{e, k^) 
is the area under the solid curves n{0p) in Fig. [51 taken 
in the angle domains < Op < n and tt < Op < 27r. 
Evaluating the integral we obtain: 



Sie,k,)^-n—± 



2A:_Ld + sgn(/i)p^ln 




where the interlevel spacing is 



ujs = uj [ I + sgn{h) 



1 p2 

irOth k±d 



(75) 
(76) 



The threshold angle 0th is of the order of £_/{2d). 

With increasing of the distance from the vortex to 
the surface the splitting of the isoenergetic lines tends 
to zero. According to the arguments presented in Sec. [TTl 
the probability of the tunneling between different qua- 
siclassical orbits can be estimated as follows [see Eq. 
dni)]: W - exp(-27r((56'p/A6'p)2), where 50p = p/^/\p\ 
and /S.6p ^ \/ y/\p\. When the splitting is so small that 
W ^ 1 (if p ^ 1), the above consideration becomes insuf- 
ficient and we should take account of the Landau-Zehner 
transitions between the quasiclassical orbits. 

To study this limit we employ the approach developed 
in the Sec. IV Bl The wave functions should vanish at the 
sample surface: 



2ir 



Jk^R(0)coa{0p-e) 



^-''>^{R{0) cos{0p-0),0p) dOp^O. (77) 



Using the same technique as in Sec. IV BJ we obtain: 



^^§-0-^]Ci{O) - r 



d e 

'""do-h 






\C2{0) = e~*''«'+2'^^'*^Ci(0), (78) 



where Cy{0) ^ C{0) and C2{0) = iD{0/h). The bound- 
ary conditions for Ci_2(^) are: 



«C2(/i7r/2) = Ci(-7r/2). 
jCi(7r/2) = C2(-W2)- 



(79) 



It is important to notice that Eqs. (jTS]) are valid until 
k^d\h\ ^ 1. Since the case ft, = corresponds to the 
vortex positioned at the focal point of a concave surface, 
the above condition means that the vortex distance from 
the focal point is much larger than the atomic scale. In 
case when the vortex is situated at the center of the sur- 
face curvature [d = — 2F, h = 1 and p = 0) the Eqs. 
(|78|) appear to be very similar to the Eq. (15) in Ref. [ij 
obtained for a vortex at the center of a superconducting 
disc. The only difference is caused by the absence of the 



15 



quasiparticlc scattering at the opposite end of the trajec- 
tory since we assume the half-infinite superconducting 
sample. Hereafter considering the concave surface we 
focus on the case when the vortex is shifted from the 
curvature center towards the boundary at the distance 
exceeding the atomic length scale: \p\ ^ 1. Then, at 
large angles 9 ^ J/{LO\/\h\p) the rapidly oscillating right 
hand side in Eq. ((78|) can be neglected and we obtain the 
solution corresponding to the case of a single vortex in a 
bulk superconductor: 



C = (Ci,C2) = (aie^ 



ieO /u. 



,0-26 



— ieO/Luh\ 



(80) 



where ai, 02 are arbitrary constants. In order to ob- 
tain the transfer matrix X in the equation C{9th) = 
XC{—0th) we need to consider the domain of small an- 
gles. Introducing new functions 



Bi{e) ^C'^(6»)e»fc-Lrfe~*'''''/2-ieie/" 
B^ie) ^ -sgii{p)C2{e)e 



ikj_d ipe-' /2~-iEie/iAj 



where ei = e{h + l)/(2/i), and a new coordinate 9 = 
\f--p{9 -I- 9q) with 6*0 = s{h — l)/(2p/ia;), we obtain the 
system of equations for Bi{9), B2{9). This system coin- 
cides with the Eqs. ([M]) . though for h < the coordinate 
X is imaginary and Eqs. (|64p should be solved along the 
imaginary axis. Using the solution of Eqs. ((64)) (see Ap- 
pendix |B]), we obtain the transfer matrix for h > 0: 



X = e-^f'^'^i - i (o-^Imn + ayRcTi) , 
where / is the unity matrix, 

n = V2sinh(7rpV2)e-"P'/^e'>'i , 

*pV2)) 



(81) 



and xi = 2kA_d + p^lii\9*\ + arg(r(l 
Analogously, for /i < we get: 



X = e"P'/2 J _ (ayReT2 + <7jmT2) 



7r/4. 
(82) 



where 



T2 = v/2sinh(7rpV2)e 



■^P^/4:„iX2 



andx2 = 2fc_Ld-p ln|6'*|-arg(r(l-ipV2)) + 7r/4. We 
denote here 9* = ^/\p\{9th + 6*0) ■ 

The quasiparticlc spectrum is obtained by matching 
the solutions in different angular domains using the trans- 
fer matrices (|8T|) . (|82p and imposing the boundary con- 
ditions dUl): 



cos 




where 



X = 2kj^d + sgn(/i) 



p^\n 



\pm 



th + VOI 



+ arg(r(l-^ 



(84) 



The energy e enters both the left- and right- hand sides 
[via 9q = e[h - l)/[2phijj)] of the Eq. ^. For small en- 
ergy (|e| ^ Ao) we have 9th ^ 9^ and the logarithm in 
Eq. ((84)) can be expanded as follows: Ind^t^-l-^olv^) ~ 
\\i{9th\/\p\) — 9q / 9th. The Landau-Zehner transitions be- 
tween quasiclassical orbits are important for i? < 1 and, 
therefore, in this limit the energy-dependent term in x 
can be always neglected: p^9o/9th 'C 1. In the oppo- 
site limit the probability of Landau-Zehner transitions 
vanishes and Eq. (|83p yields the spectrum in the form: 



2n± 



2k^d + sgii{h)p^ ln(9th^/\p\/F 



ff^Z 



(85) 



Here n is integer and wc have used the asymptotic for- 
mula for the argument of gamma function at p ^ 1: 
arg(r(l - ip^/2)) w -(pV2)[ln(pV2) - 1]. This result 
coincides with the spectrum obtained from the Bohr- 
Sommerfeld quantization rule S{e, k^) = 27rn, where 
S{e,kz) is given by Eq. ((75)) . The general scenario of 
spectrum transformation given by the Eq. ()83p is shown 
in Fig. M 
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FIG. 9: The quasiparticle spectrum for the vortex near the 
flat surface for d = 2.5^ (a), d = 1.75^ (b). The CdGM 
spectrum is shown by the dash lines. The vortex core profile 
for a single vortex is approximated by Eq. (|40|) with ^v = ^, 
kp^ = 200. 

Comparing Figs. [5] and [7] one can see that the spec- 
trum of a vortex approaching the surface is analogous 
to the spectrum of a two-vortex system where the part 
of energy branches [corresponding to the upper or lower 
sign in Eq. ()68p ] is omitted. To clarify this result we 
note that the spectrum transformation in these two sit- 
uations is bringing about by the coupling of trajectories 
with opposite momentum directions, although the origin 
of coupling is different. It is caused by the normal re- 
flection from the boundary when vortex is situated near 
the sample surface. For two-vortex or vortex-antivortex 
systems the coupling occurs due to the intcrvortex tun- 
neling and precession of trajectories around the different 
vortex cores. As we have noticed above in the special 
case of a vortex near a flat boundary the spectrum ex- 
actly coincides with the spectrum of vortex-antivortex 
system if we omit the energy levels corresponding to the 
even wave functions satisfying d'^{Q,y)/dx — 0. 
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VII. DENSITY OF STATES 

Using the results of the analysis of quasiparticle spec- 
trum transformation under the magnetic field change, we 
study the corresponding DOS modification. The calcula- 
tion of DOS can be done analytically in the quasiclassical 
limit, i.e., when the energy can be written as a function 
of the classical impact parameter b = —iJ,/k±, the trajec- 
tory orientation angle 9p and the momentum projection 
kz'. s = e(fi/k±,9p,kz). Solving this equation for b we 
find a set of isoenergetic curves ^i{9p,e,kz). Hereafter 
we focus on the consideration of DOS contribution com- 
ing from a single isoenergetic curve implying that the 
summation over the index i enumerating different curves 
should be done. 

For each isoenergetic curve the set of energy bands 
£n{kz) can be found from the Bohr-Sommerfeld quanti- 
zation rule: 

S{£n,kz)^2Tm, (86) 

where S{£,kz) is the area under the curve ^(9p,e,kz). 
The density of states (per unit length and per spin pro- 
jection) is then given by a standard expression: 

-t-fcF 



E 



dsnikz = qnis)) 



dk. 



(87) 



where the energy spectra e„(fcz) are even functions of 
momentum due to the symmetry of the BdG equations 
with respect to the z-axis inversion and (7„(e) is a set 
of positive momenta satisfying the equation Sniln) = £• 
We also assume here that \S{e,kz)\ is a monotonic func- 
tion oi kz > reaching the maximal value at kz = 0. 
This condition guarantees that we get a single positive 
Qn root for each energy branch. Such assumption appears 
to be justified for particular spectrum examples consid- 
ered above. 

Considering the differential of the function S{e, kz) for 
a fixed n index we find a simple identity 

denikz) dS/dkz 



dkz 



dS/de 



which allows us to evaluate the derivatives in Eq. ([87|l . 
As a first step we neglect the discreteness of the energy 
spectrum and replace the sum over n in Eq. ([87]) by the 
corresponding integral. Taking a fixed energy in Eq. (|86p 
one can transform the differential dn as follows: 

an = t^tt;— afcz ■ 
27r okz 

Finally the expression for DOS reads: 



v{e) = 



1 

4^ 



kp 



— kp 



dS 



de 



dk. 



Taking the spectrum of a singly-quantized vortex as an 
example we put /i = —e/uj and obtain: |S'(£)| = 27re/a;, 



v{e) = vq 



1 

2^ 



dkz 



kp 



(89) 



For a doubly-quantized vortex with the spectrum ([T]) and 
/ii_2 = —e/oj ± ^* we get v{e) = kp /['^.ujq) = 2vq. 

Now we proceed with the calculation of DOS for a two- 
vortex system and vortex near the boundary. For the 
two-vortex system the quasiclassical expression for the 
area under isoenergetic lines has the form (P5|) . There- 
fore, \dSi_2lde\ = 27r/[j and does not depend on the dis- 
tance between vortices. Thus, the DOS of the two-vortex 
system is a conserved quantity which docs not depend on 
the intervortex distance: i^(e) = kp/2u!o = 2uq. 

Considering the DOS for a vortex near the boundary of 
superconductor we use Eq. ([75]) and obtain: \dSi^2lde\ = 
27: /ojs- The important point is that lOs depends now on 
the distance from the vortex to the boundary and on the 
characteristics of the surface (i.e., the focal distance F). 
The resulting low-energy DOS can be written as follows: 
v = Vf) + sgi\{h)5v, where 



8v/v(i 



2llIq 



'K'^Qthkpd J_^.^ ujkj_ 



fcF p2 



dkz , 



(90) 



where 



Ao 



-Ko{d) 



UJ K^k^d\l + d/2F\ 

The ratio 5v /vq is a monotonically decreasing function 
of the vortex distance to the surface d and at d » ^ it 
can be evaluated as follows: Sv/vq ~ (^/d) exp(— 4(i/^). 
Therefore, if /i < (/i > 0), then the DOS is reduced 
(increased) as the vortex approaches the surface. If the 
vortex is very close to the surface {d < \F\), then h is 
always negative, therefore the DOS is suppressed. For 
example, at d = ^ the correction of DOS 61/ is of the 
order of O.li^o- This effect can be interpreted as the dis- 
appearance of the anomalous spectrum branch occurring 
as the vortex approaches the boundary and finally leaves 
the sample. The DOS reduction is the direct consequence 
of the increase of interlevel distance cus ([75)1 in the vortex 
spectrum due to the appearance of the prohibited do- 
main of trajectory orientation angles shown in Fig. HJa). 
The decrease of the distance d should result in shrink- 
ing of the quasiclassical orbits in {fi, 9p)-spa,cc and the 
DOS suppression till to the complete disappearance of 
the anomalous spectrum branch at the moment of vor- 
tex exit. 

If the spectrum discreteness can not be neglected, the 
above quasiclassical calculation of DOS becomes insuffi- 
cient. In this case, the rigorous calculation of DOS on 
the basis of the expression for the quantized spectrum 
of the two-vortex system (p5|) and the vortex near the 
boundary (j83p can be done numerically. The results are 
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shown in the Figs. [10] and [TT] To avoid the singulari- 
ties the DOS is averaged over the small energy interval 
O.lwo- In real experimental conditions such smearing of 
DOS can be caused, e.g., by finite temperature or scat- 
tering effects. 
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FIG. 10: The density of states for a two-vortex system for 
a = 5^ (a), a = 3.5^ (b). The doubled CdGM DOS is shown 
by the dash lines. The vortex core profile for a single vortex 
is approximated by Eq. (|40p with ^v = (,, kpS, = 200. 
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FIG. 11: The density of states for a vortex near the flat surface 
for d = 2.5C (a), d = 1.75^ (b). The CdGM DOS is shown by 
the dash lines. The vortex core profile for a single vortex is 
approximated by Eq. (|HJ)) with ^„ = ^, kpS, — 200. 



The DOS of the two-vortex system (Fig. \TU^ in the 
case a < Qc consists of two sets of small peaks shifted 
by the value wo(2x — [2%]) with background level of 2i/o 
[Fig. fTUTb)]. Here square brackets denote the integer 
part and x is given by the Eq. ((^ . These peaks are van 
Hove singularities corresponding to the extrema of the 
spectrum branches at fcz = [Fig. [Zl^b)]. As the distance 
between vortices increases the DOS tends to the doubled 
value of CdGM DOS of the isolated vortex, shown by the 
dash lines [Fig. [TUi;a)]. 

The expression for the spectrum Eq. (|83p of a vor- 
tex near the boundary is analogous to the spectrum of 
the two-vortex system, where part of the branches, cor- 
responding to the upper or lower sign in Eq. ([55)1 is 
omitted. The function i'{e) for the vortex near the flat 
boundary is shown in Fig. 1111 One can see that similarly 
to the case of a vortex pair the DOS reveals two sets of 
peaks, but the energy scale of DOS oscillations is now 



larger than ujq. 



VIII. HEAT CONDUCTANCE 

In this section wc calculate the heat conductance of 
vortex states in a mesoscopic superconductor focusing 
on the low temperature limit T <^ Ao when the trans- 
port is dominated by the contribution of subgap levels. 
We consider the ballistic regime and neglect the scat- 
tering effects on the boundaries between superconductor 
and normal metal leads. The expression for the thermal 
conductance reads^^: 



1 



AirhT^ '^ 



cosh^(e„/2r) 



de^ 



dK 



dk. 



(91) 



Let us introduce the function N^e) giving the number of 
energy branches crossing a certain energy level e: 



N{e) 



Y. fdkJie 

n { 



rn(fc.)) 



de„ 



dk. 



(92) 



Then the expression ([5T|) can be rewritten as follows: 



AirhT^ Jo cosh\e/2T) 



de . 



(93) 



In the temperature interval Wq < T <C Ag the dis- 
creteness of the spectrum can be neglected. In order to 
evaluate the number of states N{e) in Eq. ((93)) we use 
the quasiclassical theory, assuming that the probability 
of Landau-Zehncr tunneling between different quasiclas- 
sical orbits is small. Generally, the quasiclassical theory 
is valid only within the momentum interval k^ < k*, oth- 
erwise the interband Landau-Zcner transitions can not 
be neglected. As a result, in Eq. ((92|) we should take the 
upper limit of integration k* instead of kp. The value 
of the threshold momentum k'^ can be estimated from 
the condition that the Landau-Zehner tunneling proba- 
bility W [see Eq. ([7])] is equal to a certain threshold value 
Wth < 1- Using the Bohr-Sommerfeld quantization rule 
(I8H1) we find: 



Nis) 



k* 



' dS{e,kz) dkz 



dk. 



27r 



|5(£,o)-g(£,fc;) 

27r 



(94) 

To evaluate the integral ([93]) we consider the Taylor ex- 
pansion |S'(e,0) - S{e, fc*)| = Y.'^^o S'(")e"/n!, where 



5'(") 



de" 



|^(e,0)-5(e,fc:)| 



e=0 



As a result, we find the expansion of the heat conductance 
in power scries of T: 



T 
4:TTh 



Y^ ^n5'(")T" , 



(95) 



n=0 
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where An are expressed in terms of Rihmann function 

C(n): 



An — 



2(7i + 2)(n+l) 



1-2 



-(n+l) 



C(n + 2). 



Consequently the effective number of conducting modes 

Ny = k/kq is: 



^^ = irjE^"^^"^^"' 



(96) 



For a singly-quantized vortex we get N{e) ~ s/luq and 

27C(3) T 



N„ 



27r^ luq 



which coincides with the expression obtained in Ref. \l4 . 
For a doubly-quantized vortex N(e) = 2/x*, where fi* ~ 
kp^ and Ny ~ fi* . Now we proceed with calculation of 
the heat conductance for the two-vortex system and for 
the vortex near the sample surface. 



A. Two— vortex system 

The results of our numerical calculation of the number 
of conducting modes as a function of temperature on the 
basis of Eq. (|9ip with the spectrum ((68)) are shown in 
Fig. [T2] . The suppression of Ny a.t T -^ cjq is caused by 



120 




FIG. 12: Temperature dependence of the number of conduct- 
ing modes N^j for a two-vortex system. Curves are plotted 
for a = 2^ to a — 5^ with the step 0.5^ (from top to bottom). 
The vortex core profile for a single vortex is approximated by 
Eq. (gD]) with ^„ = £_, kpS. = 200. 

the minigap in the spectrum. One can see that at T > ujq 



the function Ny{T) grows linearly with T. Extrapolating 
this linear dependence to T = we find the residual 
number of modes TVo, which is plotted by the solid curve 
in Fig. [T3]as a function of the intervortex distance a. 



160r:r 



120 




4 a/^ 6 



FIG. 13: Residual number of modes as a function of the in- 
tervortex distance. Solid line shows the result of the exact 
calculation based on Eq. (|9ip . while dash line is obtained 
from the analytical approximate expressions (|96|l and (|47p . 



The quasiclassical procedure suggested in the begin- 
ning of this section allows to obtain a good analytical 
approximation describing the behavior of the residual 
number of modes for a < ac- To get such approxi- 
mation we evaluate the temperature independent term 
A^o = N{0)/2 in expansion ([96]) which dominates for 
rather small intervortex distances a < ac- We use the 
expression ([M)) for the number of states N{e) where the 
area confined by the closed orbits in (p, 9p)-spa,ce is given 
by (j47p . The resulting number of conducting modes vs in- 
tervortex distance is shown in Fig. [13] by the dash curve. 
Here we choose the threshold probability Wth = 0.54 
to obtain a reasonable fit to the numerical results (solid 
curve) at a ~ 2^. The critical distance ac is defined by 
No{ac) — 0. The chosen value of the threshold prob- 
ability corresponds to the critical intervortex distance 
ac ~ 4.5^. 

For well-separated vortices, i.e. when a > 2^ the split- 
ting of energy branches is small and one can use an ap- 
proximate expression (|7ip for the area S{e, k^) to obtain: 



2 

A^o = - [kp - fcl) + — In (eth^kpa/pA 



where 



VWU) 



VK/e.)2 + 4-2 



(97) 



(98) 



This expression coincides with the estimate ([T0|) in the 
limit a 3> Ct)- 
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B. Vortex near the sample boundary 

The calculation of the number of conducting modes 
for the vortex near the boundary can be carried out sim- 
ilarly to the above analysis of the two-vortex system. 
We restrict ourselves to the case when the vortex is situ- 
ated not very close to the boundary: d > ^, when we can 
neglect the distortion of vortex core profile due to bound- 
ary effects. Taking for example a flat surface {h = —1, 
F = oo) we use the Eqs. (|91|) . ((83l) and calculate the 
residual number of conducting modes Nq obtained from 
the extrapolation of the linear parts of N^ (T) dependen- 
cies to T = 0. In Fig. [T3]we plot the resulting dependence 
No{d) (solid curve). 



distances. Namely, the residual number of conducting 
modes No{d) = Ny{d,^ 0) is defined by the number of 
states at the zero energy level N{0) [see Eq. ([M]) ] which 
is proportional to the area enclosed by the quasiclassical 
orbits in (fi, 9p)-sj>ace. Therefore the shrinking of closed 
quasiclassical orbits at d '^ ^ results in the decrease in the 
NQ{d) value. At small distances to the surface d < ^t, our 
approach does not work, but it is natural to expect fur- 
ther suppression of the subgap conducting modes number 
down to zero which accompanies the vortex exit from the 
sample. 




FIG. 14: Residual number of modes A*'o as a function of the 
distance to the surface d. Solid line shows the result of the ex- 
act calculation based on Eq. (|91|l . while dash line is obtained 
from the analytical formula (|99|l . The vortex core profile for 
a single vortex is approximated by Eq. (|40p with £,v ~ ^, 
kp^ = 200. 

The estimate of the residual number of conducting 
modes within the quasiclassical theory yields the follow- 
ing result: 



IT 

where 



~2 

-sgn{n)^\n{etu^\p\/p) , (99) 



^MW 



1-1 



VWW 



1-1 



(100) 



for the particular vortex core model (|40p . Taking dc — 
2.25.^ we plot this approximate expression ([M)l for No{d) 
in Fig. [14] (dash curve) . 

One can see that the normal reflection at the sur- 
face leads to the essential increase in the residual num- 
ber of conducting modes with the decrease in the dis- 
tance d. This effect is a consequence of the minigap 
suppression. The nonmonotonic behavior of iVo(d) at 
the distances d ~ ^ has the same origin as the de- 
crease of the zero-energy DOS, which occurs at these 



Magnetic field dependence of thermal 
conductance 



To illustrate our analysis of the heat transport in the 
vortex state it is useful to consider the magnetic field de- 
pendence of the heat conductance caused by the transfor- 
mation of the vortex structure. For this purpose we plot 
schematic dependencies of the sample magnetization and 
effective number of modes at a certain finite temperature 
T > Wo in Fig. [T5l Different branches of the magneti- 
zation curve shown in Fig. [TST a) correspond to the dif- 
ferent number of vortices in the superconducting sample. 
Starting from small magnetic fields H <C Hc2 the sample 
is in the Meissner state, i.e. the number of vortices is 
zero. The number of conducting modes Ny = k/kq [see 
Fig llSf b)] is determined by quasiparticle states with en- 
ergies above the superconducting gap Aq and therefore 
is exponentially small Ny = TV^") ~ {kpL)'^e~^°/'^ pro- 
vided T ^ Tc, where L is a characteristic transverse size 
of the sample. In increasing magnetic field the supercon- 
ducting gap is suppressed leading to a slightly growing 
Ny. When the magnetic field becomes large enough to 
introduce a vortex into the sample the number of con- 
ducting modes jumps to the value N^^^ ~ T/ujq simul- 
taneously with the vortex entry. This increase in Ny is 
caused by the appearance of subgap quasiparticle states 
localized within the vortex coroi^. The next jump in 
the number of conducting modes occurs together with 
the second vortex entry. If the sample geometry fa- 
vors the formation of a giant doubly-quantized vortex, 
the number of conducting modes rises up to the value 
N 



(2) 
MQ 



(kp^). At the interval of magnetic fields where 
the giant vortex is stable, Ny is almost constant. The 
decay of the giant vortex into singly-quantized vortices 
in decreasing magnetic field is accompanied by the de- 
crease in Ny up to the value NJ^l^ > 2N^^'> (see Section 
|V]for details). While the distance between the vortices 
grows, they approach the sample surface and the number 
of conducting modes increases again, as it was shown in 
Section IVII This increase is cut off at a certain value 



N. 



(2) 



at the field corresponding to the vortex exit. 
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FIG. 15: Schematic plots of magnetization (a) and effective 
number of modes contributing to the heat conductance (b) 
in a mesoscopic sample vs applied magnetic field. Different 
branches correspond to the states with different number of 
vortices trapped in the sample. 



IX. SUMMARY 

To sumniarize, we suggest a description of a subgap 
quasiparticle spectrum in the multi-vortex state of a 
mesoscopic superconductor. Considering multiquantum 
(giant) vortices we have obtained a general analytical ex- 
pression for the quasiparticle spectrum, which is valid for 
any value of vorticity and arbitrary vortex core model. 
Taking the simplest example of the doubly-quantized 
vortex we have considered the evolution of the anoma- 
lous spectrum branches which accompanies the splitting 
of the doubly-quantized vortex. Considering the limit 
of well-separated vortices we have found the spectrum 
of vortex clusters bonded by the quasiparticle tunnel- 
ing and have investigated the crossover to the Caroli-de 
Gennes-Matricon spectrum of isolated vortices. We have 
shown that the minigap in the quasiparticle spectrum is 
absent for the intervortex distances a < a^"^ ^ln(fci?^). 
In mesoscopic superconductors it is necessary also to take 



account of the normal reflection of quasiparticles at the 
sample surface. We have shown that the spectrum of a 
single vortex placed near the parabolic surface is trans- 
formed analogously to the two-vortex system and the 
minigap in the spectrum is suppressed when the dis- 
tance from the vortex to the surface is less than the 
critical value: d < dc ~ (5/2) ln(fci?^). When the dis- 
tance is of the order of the vortex core size, the inter- 
level spacing in the vortex spectrum becomes larger than 
the CdGM value. This effect leads to the disappearance 
of the anomalous spectrum branch when the vortex ap- 
proaches the surface. 

We have analyzed the quasiparticle density of states 
and the heat conductance along the magnetic field, which 
are determined by the anomalous branches of the quasi- 
particle spectrum. At the temperatures T 3> wq ne- 
glecting the discreteness of spectrum we have obtained 
a general expression for the DOS and heat conductance 
through the characteristics of the quasiclassical orbits in 
(/i, 0p)-space. Applying the general formula to the vor- 
tex pair we have observed a significant decrease of the 
heat conductance as a function of the growing intervor- 
tex distance. Even in the limit of the zero intervortex 
distance, i.e. for a doubly-quantized vortex, the num- 
ber of conducting modes -/V„ ~ /c^^ appears to be much 
less than the value (/cfO^i which determines the num- 
ber of conducting modes for a normal metal wire of the 
radius ^. At nonzero intervortex distances and in the 
temperature region wo < T ^ Aq the effective number 
of transport modes is a linear function of temperature: 
N^ ~ Nq + aT. The splitting of the doubly-quantized 
vortex is accompanied by the decrease of the residual 
number of modes Nq and at rather large intervortex dis- 
tances we get the doubled heat conductance of a single 
vortex: Ny '^ T/ujq <^ kp^. 

Also we have shown that the normal reflection at the 
surface of the sample leads to a considerable increase in 
the heat conductance along the magnetic field when the 
distance from the vortex to the sample boundary be- 
comes rather small: £, '^ d < dc- The exit of a vortex 
from the sample is accompanied by the disappearance of 
the anomalous spectrum branch and, therefore, both the 
heat conductance and the DOS arc suppressed aX d < £,. 
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APPENDIX A: DERIVATION OF EQ. ((59)) 

At first let us prove the following formula: 



^-ik^ /2+ikx 



F{k) dk = V2^i 



^t{x^yf/2 



I{y)dy, 
(Al) 



APPENDIX B: TRANSFER MATRICES 

Let us eonsider the following system of equations: 

d 
IttBi +xBi ^pB2 , 
ox 



d 
IT—B2 - XB2 = pBi , 
ox 



(Bl) 



where ffx) is a smooth enough function defined at —00 < , -^n^- a- 4- ^ ^ ■ 

^ ' , ^,, . „ roc ,1,^ „, N , T , , , where n > and X IS a coordmate along real or imagmary 

cc < 00 and F[k) = 2-k }_^e "'^/(xjdx. Indeed, the 



integral in the right hand side of Eq. (jAip can be written 
as follows: 



J(x-yY/2 



f{y) dy 



J{x-yY/2+i. 



''yF{k)dkdy. 



axis. We start with the case of real x. Solutions of Eqs. 
(|B1[) can be expressed in terms of the parabolic cylinder 
functions D (Ref. [TqI ) with arbitrary constants di and 
da: 

Bl = diD,p2/2 xJ- +d2D,p2/2 -xJ- , 



Noting that {x - y)'^/2 + ky = {y + k - xf 12 + kx - 
fc^/2 we integrate over the y variable using the formula 



B2 



f: 



dx = ^/^^i and obtain 



P 
^/2l 

doD. 



diD. 



ip^/2-1 X 



(B2) 



e'(^-y)'/-^+i''yF{k)dkdy 



2^.pV2-i I -^\l - 



/27ri 



^-ik^ /2+ikx 



F{k)dk. 



The asymptotic expressions for the obtained solutions for 
X 3> max(l,p) are following: 



D 



Taking Eq. ^^ in the form 



V/2 \X\I- 



^sp 



-2iip 



^ /k±a f 



i^aA Xe'i'^'Cf^dfi 



^-if. /fc_L^ , g 



i^oi^/^) +ija,] X*e~''"''c*dfi 



(A2) 



D 



ipy2 



12 

-x\l - 



'2-K- 



Jx'^ /2+i{p'^ /2) ln(V2a;)+7rpV8 



Ap2/2-l ^ 



^ia;^/2+i(p^/2) ln( V22:)-3irp^/8 



3-ix^/2-i(p^/2) ln(\/2a;)-irp^/8 



r(l - ip^/2) 



we multiply it by e^^^°'^^ ^^^ 1'^ and integrate over Q^. y V i 

Using Eq. (jAip we can transform the above integrals as 

follows' where F is the gamma function. Then, we find the scat- 
tering matrix X\ coupling the solutions B = (i?i, B2) at 

» a; > and a; < 0: 

gifeia(e-ei)V4g»MVfciag»A«ei^ j^^ ^ ... 

'^ B(a; > 0) = XxBix < 0) 

A'<^±"-/'i)[is-eif-{ei-e2f](j/0\^Q^Q^ ^ in the following form: 

k±a where / is the unity matrix. 



(B3) 



Jkj_a[{0-eif/4-e''/2]-itj.''/k^a^iliei 



c^d^i 



Ti = v/2sinh(7rp2/2)e' 



-7rp^/4 ixi 



^27ri 



^^k^a[(e-elf /i+(ei-e2Y /i-Bl/2](jr0^^Q^^Q^ ^ and xi = x^ -^v^\w\^x\ +argr(l -p2/2) + 7r/4. 

Next, we consider the case of imaginary coordinate x. 
Introducing a new variable y = ~ix wc obtain an ana- 



v/2lA 
k I a 



^ik^ae /'^c*{-9) lytical continuation of the solutions ((B2)) : 



Making use of these expressions the derivation of Eq. 
([55)1 from Eq. (jA2[) is straightforward. 



Bl = diD,p2/2{yV2i) + d2Dip2/2{-yV2i) , 



Bo 



diDip2/2_i{y^2i) - d2ApV2-i(~y v 2i) 



The asymptotic expressions for this solutions at y ^ 
niax(l,p) have the form: 

Dip2/2-i{yV2i) ~ 0, 

_^iy^ I2^i{j? 12) ln(V2./)+V/8 



ApV2-i(-2/v2i) 



r(i-V/2) 
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and the transfer matrix is 

I2 = e"P'/2j + (aj^Rcra + (jjmrz) , (B4) 

where 



T2 = V2sinh(7rp2/2)e''P /'^e*'^^ , 
and X2 = y^ - V^ In |\/2j/| - argr(l - 'p^ jl) + 7r/4. 
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